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Abstract: In this pedagogically structured article, we describe a generalized harmonic 
formulation of the Einstein equations in spherical symmetry which is regular at the origin. 
The generalized harmonic approach has attracted significant attention in numerical rela- 
tivity over the past few years, especially as applied to the problem of binary inspiral and 
merger. A key issue when using the technique is the choice of the gauge source functions, 
and recent work has provided several prescriptions for gauge drivers designed to evolve 
these functions in a controlled way. We numerically investigate the parameter spaces of 
some of these drivers in the context of fully non-linear collapse of a real, massless scalar 
field, and determine nearly optimal parameter settings for specific situations. Surprisingly, 
we find that many of the drivers that perform well in 3+1 calculations that use Cartesian 
coordinates, are considerably less effective in spherical symmetry, where some of them are, 
in fact, unstable. 



1. Introduction 



Solving Einstein equations numerically is a notoriously difficult task. After many years 
of research, several well-posed formulations of the Einstein equations have been proposed 
and tested. These include constrained Arnowitt-Deser-Misner (ADM) [Q, ^, hyperbolic 
Baumgarte-Shapiro-Shibata-Nakamura (BSSN) Q and characteristic evolution Q, just to 
name a few: we refer the reader to |6|, Q for reviews of these and other approaches. 
Among the ingredients that are key to the success of any particular formulation are 1) 
an appropriate choice of dynamic variables that results in a well-posed system, and 2) 
a choice of coordinates that remain regular during the course of the evolution. In this 
paper we focus on a specific well-posed approach known as the generalized harmonic (GH) 
formulation. This form of the Einstein equations has recently attracted significant attention 
in the numerical relativity community, in large part because of its use in obtaining the first 
long-term evolution of binary black-hole inspiral and merger |^ |9|, |l^ . 

In essence, the GH approach is a way to write the field equations such that the resulting 
system is manifestly hyperbolic, taking the form of a set of quasi-linear wave equations for 
the metric components. The basic idea underlying the strategy has a long and distinguished 
history: specifically, the use of harmonic coordinates has been instrumental in establishing 
many fundamental results in General Relativity (GR) including the characteristic structure 
of the theory |11], and the well-posedeness of the Cauchy problem for Einstein's equations 



12, |T^. However, from the computational point of view, harmonic gauge^ can be too 



restrictive, and numerical implementations using it may develop coordinate pathologies, as 



described, for instance, in and [15|. More recently, it was realized by Friedrich [16|, 



and independently by Garfinkle [17|, that much of the coordinate freedom apparently lost 
by the specific choice of harmonic gauge could be regained through the introduction of 
certain gauge source functions, while at the same time maintaining the desirable property 
of strong hyperbolicity of the field equations. In fact, the source functions can be thought 
of as representing the coordinate freedom of the Einstein equations, and when construct- 
ing solutions of the equations, via an initial value approach, for example, they must be 
completely specified in some fashion. 

Following Garfinkle's pioneering use of the generalized harmonic approach in his study 
of generic singularity formulation in cosmologies with scalar field matter , the technique 
was successfully employed by Pretorius |9|, , and subsequently by others [|^, |l^, ^ , 
for simulations of binary black hole coalescence. However, the total number of physical 
scenarios studied so far using the GH approach is limited, and there is an argument to 
be made for a more systematic exploration of the method's potential. This is especially 
the case given the relative lack of proven prescriptions for choosing the gauge functions 
appropriately in instances where the gravitational field is highly nonlinear and dynamic. 
Moreover, in order to expedite experimentation with the approach, we feel that it is useful 
to start with systems with a high degree of symmetry. Restriction to highly symmetric 
spacetimes reduces the effective spatial dimensionality of the partial differential equations 
that must be solved, yields algebraically simpler equations, and, overall, leads to enormous 



^In this paper "gauge" means "coordinate choice", and we use both expressions interchangeably. 
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savings in the computational resources required to simulate a single spacetime. This in turn 
allows for much more detailed and thorough surveys of the multi-dimensional parameter 
spaces that typically arise from a given choice of gauge functions. 

In this paper, then, we focus on the application of the generalized harmonic approach 
to the problem of gravitational collapse in spherically symmetric D-dimensional spacetime. 
Even with the restriction to spherical symmetry, we find that the strong-field aspects of the 
collapse process present significant challenges regarding the choice of the gauge functions. 
Ironically, some of these challenges may in fact be related to the symmetry restriction 
itself. As usual, in situations where a black hole forms, care must be taken to avoid the 
central singularity. This can be done through the use of singularity- avoiding coordinates, 
by excising the singularity from the computational domain, or with a combination of both 
strategies. Within the context of the GH formulation any such strategy must also be 
preserve the strong hyperbolicity of the field equations. 

Although we view our study of the GH approach for spherically symmetric collapse 
as interesting in its own right, a primary goal of this research is to prepare for an inves- 
tigation of axially symmetric systems using an analogous formulation. We thus consider 
our spherically-symmetric set up as a valuable toy model for the phenomenologically richer 
axisymmetric situation. In both cases it is natural to use coordinates in which the sym- 
metries of the spacetime are explicit. These coordinates, however, are formally singular: 
at the origin in spherical symmetry, and on the axis in axial symmetry. Thus, in both 
instances the field equations have to be regularized in numerical implementations, and one 
of the results of our work is a regularization procedure that is compatible with the GH 
approach. Moreover, we expect that the experience gained from our spherically symmetric 
calculations concerning how to choose gauge source functions will also prove useful for the 
more general case of axisymmetric computations. 

In order to maximize the usefulness of this paper to other researchers interested in 
experimenting with the generalized harmonic approach, we have attempted to make the 
following presentation reasonably self-contained and pedagogical in nature. We thus begin 
in Sec. ^ with a brief presentation of the basic GH formulae in full generality, along with a 
discussion of the constraint equations. Although the constraints are consistently preserved 
by the GH evolution equations in the continuum limit, in numerical calculations at finite 
resolution, deviations from the constraints generically develop. In order to maintain sta- 
bility these deviations must be damped and we describe a method that effectively achieves 
this damping. Sec. ^ is devoted to a detailed discussion of coordinate conditions. One key 
issue that we consider is the non-trivial problem of prescribing the GH source functions to 
mimic some of the more popular and successful coordinate conditions that have historically 
been used in numerical relativity calculations. Following recent proposals ||9|, [l^ we 

describe the formulation of the gauge conditions as hyperbolic evolution equations: is this 
approach the gauge functions are evolved, or "driven", to desired targets in a controlled 
way, rather than being fixed instantly. 

In Sec. ^ we adapt the GH formulae to the case of asymptotically flat, spherically 
symmetric configurations in D spacetime dimensions. We derive the field equations, cast 
them into a form suitable for numerical solution, discuss initial and boundary conditions, 
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and regularize the singular origin by introducing a new variable. Since we use spherical 
coordinates adopted to the symmetry, the GH source functions appear to diverge at the 
origin as 1/r. Hence, we regularize these functions as well by subtracting off the singular 
contribution that appears in the flat spacetime limit. The operators that appear in the 
various gauge drivers then act on the regularized source functions. 

In order to endow our model with non-trivial dynamics, we introduce a minimally 
coupled, real, massless scalar field. The initial distribution of the scalar matter is freely 
specified and our results are grouped according to the "strength" of the initial data. In 
each case we simulate the time evolution of a single Gaussian pulse of scalar field that is 
initially centered at the origin. The weak and the intermediate data correspond to the 
dispersion of relatively dilute pulses, while a typical strong data configuration collapses to 
form a black hole or nearly does so. 

Mathematically, the task of treating the coupled Einstein-scalar system involves the 
solution of a set of several quasi- linear wave equations. (Here we note that some of the gauge 
drivers involve auxiliary variables that obey first-order- in-time differential equations.) Our 
numerical approach to solving this system using finite difference techniques is detailed in 
Sec. |5[ We compactify the spatial (radial) dimension into a finite region and cover it by 
a discrete lattice. This allows us to include spatial infinity on the finite difference mesh, 
which has the advantage of enabling us to set exact boundary conditions corresponding to 
asymptotic fiatness. Following we directly discretize the second-order-in time- wave- 

equations on the mesh, and use a point-wise Gauss-Seidel relaxation method to update 
the discrete unknowns at each time step. In order to damp high-frequency components of 
the numerical solution — which can generically lead to instabilities — we incorporate explicit 
dissipation of the Kreiss-Oliger type. This dissipation is also essential for attenuating 
spurious reflections from the outer region of the compactified domain that would otherwise 
quickly contaminate the solution in the interior (i.e. near the origin). 

For the case of black hole formation we have investigated both of the approaches 
mentioned above for avoiding the central physical singularity. On the one hand, we have 
implemented an excision technique, in which an excision surface is chosen so that all char- 
acteristics on it are pointing inwards, obviating the need for explicit boundary conditions 
for the evolution equations. On the other hand, we have also experimented with the use of 
singularity avoiding slicing conditions, that "freeze" the evolution in the strong curvature 
regions. However, we flnd that in our case the calculations using singularity- avoiding slic- 
ings tend not to run as long as those with excision and appear to crash prematurely due 
to numerical errors that build up in the strong curvature regions. 

Sec. ^ is devoted to a discussion of our detailed investigation of the performance of 
several coordinate conditions as applied to calculations involving various strengths of ini- 
tial data. As already mentioned, the parameter spaces associated with many of the gauge 
drivers that we consider here are multidimensional. Thus, even with the significant reduc- 
tion in needed computational resources that the restriction to spherical symmetry provides, 
we have not found it feasible to identify optimal parameters in all cases. In some instances 
then, we simply report what appears to be typical behavior for a particular gauge, while still 
trying to explore the effects of the variation of key parameters on the quality of the solutions. 
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Interestingly, we find that several of the gauge drivers that have been successfully used in 



the 3 + 1 simulations of black hole collisions that use Cartesian coordinates 10, 18, 19 1 
are considerably less effective for our spherically symmetric calculations. In particular, it 
is not always possible to drive the lapse to a certain value as reported in |l^], nor is 
it always possible to enforce a desired gauge for a long time by using one of the drivers 
described in p^ . Overall, our calculations seem to be more sensitive to the specific choices 
of parameters for the drivers than the Cartesian computations, and this is an issue which 
warrants further investigation. 

Nevertheless, our results indicate that with a certain amount of parameter tuning, 
several of the gauge conditions that we investigate facilitate the simulation of many inter- 
esting scenarios. We are thus encouraged by this particular application of the generalized 
harmonic approach, and our conclusions and discussion in Sec. |^ includes an outline of 
some future extensions of the work. 



2. Generalized harmonic formulation 

We consider the Einstein equations on a D-dimensional spacetime and written in the form 
Rfiu = ^ttGnT^i, = SttGn (^fiu - _ ^ QnvT^ , (2.1) 

where Qf^i, is the metric, -R^j^ is the Ricci tensor, T^u is the energy-momentum tensor of the 
matter with trace T, and Gat is the D-dimensional Newton constant. Hereafter, we adopt 
units for which SttGn = 1- 

The Ricci tensor that appears in the left-hand-side of ( |2.1| ) contains various second 
derivatives of the metric components g^^: these second derivatives collectively constitute 
the principal part of R^u, viewed as an operator on g^^^. This principal part can be decom- 
posed into a term g'^^ dapgixv, plus mixed derivatives of the form g'^'^ dai_igju ■ Without the 
mixed derivatives, ( |2.1D would represent manifestly (and strongly) hyperbolic wave equa- 
tions for the (7^jy [|2l|. Strong hyperbolicity is a highly desirable property since mathematical 
theorems then ensure (local) existence and uniqueness of solutions at the continuum level. 
This, in turn, means that it should be possible to construct stable (convergent) numerical 
discretizations of the field equations. 

One can view the generalized harmonic (GH) formulation of general relativity as a 



particular method that eliminates the mixed second derivatives appearing in (p. 11) pq, 17 



|8|, 10, As the name suggests, the technique generalizes the harmonic approach in which 
the spacetime coordinates, x^, satisfy the harmonic coordinate condition 

□x" = 0. (2.2) 

Here we have 

□x" = -^d, iV^g'"') = -r" = -g-'^T^p, (2.3) 

where F"^ are the usual Christoffel symbols. 



- 4 - 



It was realized by Friedrich [16| and also by Garfinkle [0], that it is possible to elim- 
inate the mixed derivatives in the principal part of the Einstein equations while largely 
recovering the coordinate freedom than is lost by choosing the harmonic gauge. Instead of 



( ]2.2| ), one requires that that the coordinates satisfy 

□x" = ff", (2.4) 

where = Qa/sH^ are arbitrary "gauge source functions" ^ which are to be viewed as 
specified quantities. One then defines the GH constraint 

= - Dx", (2.5) 



which clearly must vanish provided ( p.3[ ) holds, and then modifies the Einstein equations 
as follows: 

Rfiiy - ^^(^t;!/) = T^u- (2.6) 
This last equation can be written more explicitly as 

Now, provided that the are functions of the coordinates and the metric only, but not of 



the metric derivatives — namely = Ha{x, g) — the field equations (2.7) form a manifestly 
hyperbolic system. We reemphasize that the source functions are arbitrary at this 
stage and that their specification is equivalent to choosing the coordinate system for the 
spacetime under consideration ("fixing the gauge"). Determining an effective prescription 
for the source functions is thus crucial for the efficacy of the GH approach, and several 
strategies for fixing the are discussed in the next section. 

Having prescribed the coordinates we integrate the equations forward in time. Con- 
sistency of the scheme requires that the GH constraint ( |2.5D be preserved in time. The 
contracted Binachi identities guarantee that this is indeed the case, since, using those 
identities, one can show g H that itself satisfies a wave equation, 

aC^ + R'^^C = 0. (2.8) 

Thus, assuming that the evolution is generated from an initial hypersurface on which 
= 5tC° = 0, (|2^ ) guarantees that C° = for all future (or past) times. 
Although the GH constraint is preserved at the continuum level, in numerical calcu- 
lations, where equations are discretized on a mesh with some characteristic mesh scale, h, 
the constraint cannot be expected to hold exactly. More troublingly, experience shows that 



numerical solutions of (2.7) — particularly in strong field cases, such as those involving black 
holes — can admit "constraint violating modes" , with the result that the desired continuum 
solution is not obtained in the limit /i — >• 0. Fortunately, an effective way of preventing 



^In a slight abuse of notation and terminology we will refer to both Ha and H" as "the" gauge source 
functions. 
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the development of such modes in numerical calculations exists: one adds terms to the 
field equations that are explicitly designed to damp constraint violations (see e.g. p3|). In 
our implementation we follow Pretorius by adding constraint damping terms in a 



fashion inspired by studies of the so-called 7-systems 25 1. The modified equations take 
the form 



1^ {^{pfiu) - ^9tiiy nf^Cp^ = T^y. (2. 



Here, n" is the future-directed, unit time-like vector normal to the t = const, hypersurfaces, 
which can be written as 

na = - dat, (2.10) 

and K, is an adjustable parameter that controls the damping timescale. Specifically, as 
discussed in |25], small constraint perturbations about a fixed background decay exponen- 



tially with a characteristic timescale of order k. We note that the constraint damping 
term contains only first derivatives of the metric and hence does not affect the principal 
(hyperbolic) part of the equations. 

3. Coordinate conditions 

As we have already mentioned, fixing the coordinates in the GH approach amounts to 
specifying the source functions H^. In this regard, it is instructive to examine the rela- 
tionship between the and the lapse function and shift vector that appear in the ADM, 
or space-plus-time, formulation of general relativity. We recall that in the ADM formalism 
the line element can be written as 

ds^ = -a^dt^ + jij {dx' + P'dt) {dx^ + (3^dt) , (3.1) 

where a is the lapse function, /3* is the shift vector, and 'jij is the spatial metric of the 



t = const, hypersurfaces. Using this form of the spacetime metric in (2.4) yields 



dta — jS^dka = —a {Hn + aK) , 
dt/3' - P^duP' = ai^ \a (h, + ^""-'^Tj^n"') - djo] , (3.2) 



where Hn = n'^H^ = {Ht — (3^ Hi) /a is the normal component of the source function H^, 
K is the trace of the extrinsic curvature tensor of the t = const, slices, and the '-^"^^Tj^; 
are Christoffel symbols associated with the spatial 7ij. Bearing in mind that the temporal 
component of the source function is thus determined by Ht = aHn + P^Hi, these last 
equations clearly exhibit the connection between the gauge source functions and the time 
evolution of the lapse and shift. 

In his groundbreaking application of the GH approach ||9|, |lO|, Pretorius used insight 
derived from considering this relationship between the H^ and the ADM kinematic variables 
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to devise a methodology that generates effective gauge source functions for the problem 
of binary black hole collisions. His strategy elevates the status of the Ha to independent 
dynamical variables that satisfy time-dependent partial differential equations. Crucially, 
the evolution equations for the Ha are designed so that the lapse and shift which (implicitly) 
result from the time development have certain desirable properties. For example, the 
equation for Ht is tailored in an attempt to keep the value of the lapse function of order 
unity everywhere — including near the surfaces of the black holes — during the evolution. 

One specific prescription for achieving this type of control evolves the gauge source 
functions according to 

a — an 

UHt = -6 ^ + 6 Ht^^n^", 

Hi = 0, (3.3) 

where □ is the covariant wave operator, and 0:0,^1,^2 and q are adjustable constants'^. 
Thus the temporal source function satisfies a wave equation similar to those that govern 



the metric components in the system (2.9). The first term on the right-hand-side of ( p.3| ) 
is designed to "drive" Ht to a value that results in a lapse that is approximately oq- The 
second, "frictional" term tends to confine Ht to this value. For the case of the spatial 
coordinates, Pretorius found that the simplest choice of spatially harmonic gauge — Hi = 
— was sufficient in simulations of binary black hole collisions. Importantly, the choice 



( |3.3| ) ensures that the hyperbolicity of the combined evolution system is preserved. A 
slight generalization of this technique was considered in |1S] where instead of using Hi = 0, 
the spatial components of the source functions are evolved according to 

UHi = -iA+i2H,^^n^' (3.4) 

where ^3 is an additional parameter. 

One possible problem with the specific driver approach outlined above is that the 
coordinates that result do not correspond to those produced by any of the more familiar 



coordinate conditions typically used in numerical relativity. Recently, Lindblom et al |18| 
proposed driver conditions that are crafted so that the source functions that result imply 
particular conditions on the corresponding lapse and shift. We now proceed to a review of 
this interesting and promising approach. 

We begin by observing that many traditional coordinate conditions of numerical rel- 
ativity can be written as = Fa{x, g, dg) where the are to be viewed as "effective" 
gauge source functions which could be computed, for example, were the entire spacetime in 
hand. Within the GH approach, enforcing such a condition algebraically by simply setting 
Ha = Fa will generally destroy the hyperbolicity of the system, since the H(^^.y-^ terms in 
(|2.7|) will generically give rise to mixed second derivatives of the metric. Lindblom et al 



^Sometimes it is convenient to assume that ^1 and ^2 are given functions of space and time rather than 
mere constants. For example, one might require that the gauge driver is switched on gradually in time, 
or that it be active only in certain regions, e.g. in the vicinity of a black hole, and that its effect vanish 
asymptotically, so that pure harmonic coordinates are recovered at large distances. 
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circumvent this difficulty by generalizing (|3.3| ) to 



OHo, = Qa{x,g,dg,H,dH), (3.5) 

where O is a second order hyperbolic operator and Qa is chosen so that the source functions 
evolve towards the concrete Fa = Fa{x, g, dg) that define the desired gauge. The combined 
system ( p.9D and (|3.5| ) will remain hyperbolic provided the Qa depend on at most first 
derivatives of the fields. In analogy with (3.3) the authors of [p!^] choose 



Qa = f4iHa-Fa) + 2^i2dtHa+vWa, (3.6) 

where m,fj.2 and rj are adjustable parameters, and Wa is assumed to satisfy 

dtWa+vWa = dHa, (3.7) 

where O is the part of O that contains only spatial derivatives. When the spacetime is 
stationary, time-derivatives vanish and equations ( |3.5D and ( p. 61 ) then imply Ha = Fa- 
Notice that without the introduction of the auxiliary fields, Wa, this property could not 



be attained for general, position dependent gauges |18|. 

In order to implement this method for a specific desired gauge choice one must first 
compute the corresponding target source functions. Fa- Here we focus on gauges of the 



schematic form Ga{x, g,dg) = for which one can choose [18 



Fa = -ra-qGa, (3.8) 

where g is a tunable parameter. In the GH formalism. Ha = — r^, and ( p.8| ) then implies 
Ha — Fa = q Ga- This demonstrates that when the GH constraint is satisfied. Ha is driven 
to Fa if Ga is driven to zero. We next discuss several specific coordinate choices that are 
explored in this paper. 

3.1 Slicing conditions 

For the particular choices of the slicing conditions that we use in this paper, it is more 
convenient to calculate the normal component of the target source functions, Fn = n^F^ = 
{Ft — l3^Fi)/a, than the temporal component. Ft, itself (see (3.2)). Once this is done, then 



in conjunction with the shift conditions that fix Fi, the temporal component can be easily 
computed via Ft = aFn + 

• Constant curvature slicing, K = Kq. Here we assume that the trace, K{g,dg), of 
the extrinsic curvature of the spatial slices is constant. When Kq = 0, we have 
the famous maximal slicing condition pSj whose significant popularity in numerical 
calculations is due in large part to the strong singularity- avoiding property exhibited 
by the resulting constant-time surfaces (see Sec. The constant curvature foliation 
can be written as G„ = 0, where 

Gn = Ko- K = Ko + (3.9) 
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• Bona-Masso slicing ||30|] . This condition can be written as 



Gn = (dta - /3'dia) + a^f{a) (K - Kq) , 



(3.10) 



where /(a) is an arbitrary function of the lapse. ^ The choice /(a) = 2/a corresponds 
to the popular 1 + log slicing. 



In terms of implementing these slicing conditions, we note that 



implies 



(3.11) 



where g„ is a parameter, and that the kinematic quantities such as the lapse and shift 
which appear in various formulae above can always be written in terms of the fundamental 
dynamical variables of the scheme (i.e. the metric components and their first derivatives). 

3.2 Shift conditions 

An important class of shift conditions which is often used in numerical relativity employs 



versions of the so-called F-driver |27]. In this approach, one first introduces the conformally 
rescaled spatial metric, jij = 'j'^jij, with 7 = det7ij and a an arbitrary parameter, then 
computes the contracted Christoffel symbols. 



{D-i)f.i UD-i) Y{.;ykj 
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l + a{D-2,) 



7^^'aiog7 + 7^47^^' 



(3.12) 



and imposes certain conditions on their dynamics. The F-driver strategy is related to the 
minimal distortion condition [28, 29 1 which is designed to minimize the time variation of 
7ij (see e.g. 



F-freezing. Here one requires 

5^(^-1)^ = 0, (3.13) 
which implies that during the evolution (^~i)f* is fixed, (^-i)f* ={D-i) Fol- 



lowing [18 1 we attempt to evolve to this choice by choosing 

G, = 7iJ(^-i)P(0)-(^-i) 



F-driver. Again following [18[ we write the driver condition as 



dtB' + 5' =(^-1) r, 
where v and r] are adjustable parameters. Then one can choose 

= 7ii idt^' - y ^""^'^f ^' + vm 



(3.14) 

(3.15) 
(3.16) 

(3.17) 



^Sometimes the geometric derivative dna = {dtOL — P'' dio) / ol is replaced with the partial time derivative 



dta. 
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The auxiliary variable is evolved using ( 3.16| ) and it is important to note that 



adding this equation to the scheme does not destroy the hyperbolicity of the combined 
evolution system |18|. 

We have also experimented with a geometric version of the driver where the partial 
time derivative dt in ( |3.16D is replaced with the covariant derivative n^V^ = (dt — 

Implementation of the above shift conditions is effected by setting the corresponding 



spatial target source function defined by (3.8) according to 



Fi = -Ti-qiGi, (3.18) 

where qi is an adjustable parameter. 

4. Spherically-symmetric reduction 

Having described the basics of the GH formalism, we now specialize to spherically sym- 
metric spacetimes. We consider a D-dimensional spacetime with SO{D — 2) rotational 
symmetry, and write the Z)-dimensional line element in the form 

ds'^ = gl^^dx^dx" = g^^^dx^dx^ + e^^d^l. (4.1) 

Here is the metric on a unit n-sphere, n = D — 2, a,b = {t,r}, and the metric g^J^'' 
and scalar S are functions of t and the radial coordinate, r, alone. 

Although we will later specialize to the case of a real, massless scalar field, for generality 
we first adopt as a matter source minimally coupled complex scalar field, with a potential 
y(|$|). The action that describes the system can be written as 

S = J ^-gi^) [r^'^^ -da^d"^* -2V{\^)^ dx^. (4.2) 



By varying the action with respect to the fields one gets the Einstein equations ( |2.lD with 
the energy-momentum tensor = \ (f^^^^M*^* + "^z^*^* ^M^) + ^/(-^ ~ 2)5}^?^ V, as weh 
as the general relativistic Klein-Gordon equation for the scalar field. Specifically, the GH 
transformation of the Einstein equations as given by (^^) reads 



4S - ^(^.^.) = ^^S! ^' (4-4) 
= dVjd^*, (4.5) 

where is the D-dimensional Ricci tensor and di are the angular coordinates. In 

spherical symmetry it suffices to use any specific angular component of the Ricci tensor, 
and for convenience we use ^q^^ where Q\ is defined by dVi^ = dOf + sin^ 9idQ,'^_i. 
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The form of the metric (|4.1| ) is not yet optimal for use in numerical computations. In 
this paper we are mostly interested in asymptotically flat solutions and thus the following 
section describes a more natural ansatz for use in that instance. 



4.1 Spatial asymptotics 

In spherical coordinates, flat spacetime can be written as 

ds'^ = -df + dr"^ + r'^ dOi. (4.6) 



It follows from ( |4.1| ) that asymptotically Qab — ^ where "r^ab is a Minkowski metric, 
and S — )• log r, (i.e. S diverges at spatial infinity). Since this divergence complicates the 
numerical implementation of boundary conditions, we introduce a new function, S, defined 
by 5 = 5 — log r, which is regular everywhere. We then adopt the following, more regular 
form for the line element in the asymptotically fiat case: 

ds^ = Qabdx^dx^ + ^d^l- (4.7) 



In spherical coordinates, the source function derived from ( |2.4| ) does not vanish even in fiat 
spacetime where it becomes 



^Mmk ^ _pMmk ^ ({)^n/r, {n - 1) cot 0i, (n - 2) cot ^2, • • • , cot 0„_i, 0). (4.8) 

Since near the origin spacetime is locally fiat, the radial component of the source function 
is generically singular at r = 0, diverging as n/r. To regularize this radial component, we 
thus subtract the singular background contribution by transforming Ha — )• + JQ-ff^™*^, 
and use the functions Ht and Hr defined by 

Ha = {Ht{t,r),Hr{t,r) +n/r, (n - l)cot6li, (n - 2) cot 6^2,. . . , cot 6'n-i, 0) . (4.9) 

in our formulae. 



With the line-element (4/7) and the source functions (O), the asymptotic behavior of 
the fields is simply 

Oab^Vab, S^O, 0^0, Ht^O, Hr^O (4.10) 

In App. |A| we also analyze the asymptotically AdS spacetime, which is described in our 
model ( |4.2D for the case that the scalar field potential satisfies V{0) — t- A < 0. 



4.2 Center of symmetry, r = 



Invariance of the line element (O) under the refiection r — t- — r in spherical symmetry 
implies that gtr is an odd function of r, while gu,9rr, S and $ are even in r. Additionally, 
the GH constraint (2^) implies that the source functions Hr, regularized via (O), and Ht 



are odd and even in r, respectively. 



- 11 - 



Moreover, the requirement that the surface area of an n-sphere must vanish at the 
origin^ imphes grr{t,0) = e^'^^*'^^. We note that this is an extra condition on S, which 
thus has to satisfy both this relation, as well as the constraint that it have vanishing radial 
derivative at r = — specifically that grr — e^^ = O(r^). Therefore, at r = we essentially 
have three conditions on the two fields S and grr- In the continuum, and given regular 
initial data, the evolution equations will preserve regularity: however, in a numerical code 
that solves the equations discretized on a lattice, this will be true only up to discretization 
errors. As a general rule-of-thumb, the number of boundary conditions should be equal to 
the number of evolved variables in order to avoid regularity problems and divergences of a 
numerical implementation. 

An elegant way to deal with this regularity issue involves definition of a new variable, 

A: 6 

^_ grr-e^\ (4.11) 
r 

At the origin one then has A ~ 0(r). Therefore, after changing variables from 5 to A by 
using S = (1/2) log{grr — r X) in all equations, and imposing A(t, 0) = at the origin, one 
ends up with a system where there is no over-constraining due to the demand of regularity 
at r = 0. In addition, we note that at spatial infinity we have A = 0, and that the 
hyperbolicity of the GH system is not affected by the change of variables. 

However, as described in detail in Sec. 5^, we were able to implement a more straight- 
forward regularization method that maintains S" as a fundamental dynamical variable, and 
thus opted to use that approach in our current calculations. 



4.3 The equations 



With the metric ansatz ( ^^ ) and the regularized source function (|4.9|), equations (4.3)~ 



become 5 equations for the 5 variables, gtt, 9tr, 9rr, S and that schematically can 
be written as^ 

-lg'''9ab,cd + --- = l{da<^ db^* + da<^* 8^^) + -^^OabV, (4.12) 

9'''S,,, + --- = -^^V, (4.13) 
g'"^^,cd + --- = dV/d^*. (4.14) 

Here ellipses denote terms that may contain the metric and/or the source functions, as 
well as their first derivatives in various combinations (see App. |^ for the explicit set of 
equations in the four-dimensional case). These equations are to be evolved forward in time 
starting from the initial {t = 0) time slice, where values for the fields and their first time 
derivatives must be prescribed. 



^that is, that the radial and areal coordinates coincide at the origin, to avoid a conical singularity there. 
®We note that a similar variable was introduced in |^^, also for the purpose of regularization. 
^Using A instead of S does not change this structure since the equation that governs A is a linear 
combination of the equations that govern S and grr- 
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4.4 Coordinate choices 



Here we adapt the prescriptions for choosing the gauge functions (Ht and Hr) that were 
described in Sec. ^, to the case of spherical symmetry. We again note that the radial 
source function is singular at the origin in spherical symmetry, and that we thus regularize 
it via (^^). Since this regularization involves subtracting the flat-spacetime singular part 
from Hr, any specific coordinate conditions discussed here are thus defined relative to 
spherical Minkowski spacetime. 



For the case of the gauge condition (|3.3D inspired by Pretorius' original work, we have 
DHt = -6 + 6 {dtHt - 13 drHt) /a, (4.15) 

Hr = 0. 



Similarly for the modification of the above proposed in [ [l9| |, we have (using (3.4)) 



DHt = -Ci + 6 {dtHt - /3 drHt) /a, (4.16) 

DHr = 4 + 6 {dtHr - 13 drHr) /«. 

In the above equations □ is the regularized scalar wave operator in spherical symmetry, 
given by 

UH^ = g^^'^d^d.H^ - {t'' + g^^ '^5^^.5^) d, H^. (4.17) 

Turning now to the case of the gauge drivers introduced by Lindblom et al, we note that 
the operator in ([4. 9]) is essentially the vector d'Alambertian^ [18| 



OH^ = g'^'d^d.H^ - T^d.H^ - 2 g^^'^Vl^d^Hp + f i?^ - dj:^] Hp. (4.18) 



In order to avoid having second-derivatives of the metric, the Ricci tensor in the last term 
should be thought of as being determined by matter sources and replaced with Ta, in 
accordance with the Einstein equations. In addition, using the GH constraint Ht^ — — Tct; 
the term —daT^Hj^ is replaced with —daH^^Tp. Finally, we regularize the operator by 
subtracting the irregular contributions that appear in the fiat spacetime limit. After these 
manipulations we arrive at 

OH^ = g'^^d^d.H^- (r^ + /'^ ^6';5^) d, H^-2 g^^'^Vl^d^Hp- (f^ + do^H^) {Vp + ^5„„) 

(4.19) 

where 5^ is a Kronecker delta, and there is no summation over the index a. 

The target source function, F„, is determined by ( |3.9[ ) or ( 3.101 ), and by ( 3.11 ). The 



*-ffa does not transform as a vector undergauge transformations, so the equation should be understood 
as written in particular global coordinates ||l8|; in the current case, these are our spherical coordinates. 
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lapse and shift are given in terms of the metric components, 



a = \l -gtt + glrldrr, 
/3 = gtr/gr-r, 



(4.20) 



as is the trace of the extrinsic curvature (see ( [B.2| ) for the explicit form). 

Our shift conditions involve the contracted conformal Christoffel symbols, Fj, defined 
by ( |3.12| ), and in spherical symmetry the only non-trivial component is ^^~^^Tj. given by 



(^-i)f . = -n (1 + (n - 1) a) S' + l^^S^^U 



9r 



(4.21) 



Here ()' = dr, and we have used the fact that ^rr = grr- Once again, in order to obtain a 
regular expression we have subtracted the flat-spacetime term, {^-i)pMmk _ _j^(2^ _j_ — 
l)o")/r, which is singular at the origin. 



The target function for the F-freezing condition ( 3.14 ) takes the form 



(^-i)f,(o. 



gr 



grr{0,r) 



2na[S-S{0,r)] _{D-1) f 



(4.22) 



where = + n/r is the D-dimensional connection which has also been regularized via 
subtraction of an irregular flat-spacetime term. The explicit expression for F^ is given in 



For the case of the F-driver condition (3.18) in spherical symmetry, the target source 
function is 



gr 



(4.23) 



where an over-dot denotes partial differentiation with respect to t. The auxiliary field B 
is evolved using 

ij + r?2i?=(^-i)f,(e2«^5„,)"V5rr. (4.24) 



4.5 Initial data 

We now consider specification of initial data, which as stated previously, are values for the 
fields and their first time derivatives at t = 0. For simplicity (and without much loss of 
generality), we restrict attention to time-symmetric initial conditions. 

Given the assumption of time symmetry at t = 0, initial data for the scalar field reduces 
to the specification of <I>(0,r), which we take to have the form of a Gaussian, 

$(0,r) = ^oe-^"-"")'/^', (4.25) 

where ^>0) ''o and A are adjustable parameters. 

The momentum constraint is trivially satisfied for time-symmetric initial data, and 



-14- 



writing the initial metric as 



(4.26) 



the Hamiltonian constraint becomes a non-hnear ordinary differential equation for ^{0,r) 



n 



< + -iP' + (n - 2) 



+ 



1 /I 



2n V 2 



-^>'^>*' + V y = 0. 



This equation is solved using the boundary conditions •(/''(O, ?^)|r=o = and -0(0, r)|, 
and then once ij) has been determined, the metric components are initialized via 

9rr = i'^, 

S = 2 1ogV', 
gtr=P' = \ = 0. 



(4.27) 

X) — 1) 

(4.28) 



For time-symmetric initial data we require that all first time derivatives of the metric 
components vanish. 

We next determine the initial conditions for the lapse and the variables used in the 
gauge drivers. We begin by setting Ht{0,r) = Hr{0,r) = 0. Using 



Hr{0,r) = -r, 0,r = - + 2 n - 1^^, 



(4.29) 



we obtain an equation relating a(0, r) to the initial value of Hj.. With our choice, Hr{0, r) 
0, this equation can be integrated to yield 



a(0,r) = V'(0,r)"2("-i). 



(4.30) 



Next we require that the target coordinate conditions are initially satisfied, namely that 
Fa{0,r) = Ga{0,r) = Ha{0,r) = 0. We note that since time-symmetry implies K(0,r) = 
Kq = 0, the normal component of the gauge function for the constant curvature folia- 
tion vanishes, Gn(0,r) = QuKq = 0, as it does for the Bona-Masso slicing, G,i(0,r) = 
—Qn a(0, r)^/(a(0, r)) Kq = 0. The T-freezing condition ( 4.22 ) obviously satisfies Gj(0, r) = 
0, while requiring this for the T-driver condition ( 4.23| ) will set the initial value of the aux- 
iliary field i?^, 

5(0, r) =(^-1) f,e-2-^-5-7r?2|t=o. (4.31) 
Here the initial value for the radial component of the contracted conformal Christoffel 



symbol r,.(0,r), defined by ( 4.21 ), is found using the relations (4 



(^-i)f,(o,. 



-2(n-l)(l + (n + l)cT) 



(4.32) 



^Note th at for time-symmetric initial conditions this consistently coincides with the values of _B(0,r) 
found from (124). 
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Figure 1: The compactificd domain of integration, and the numerical lattice. Our finite difference 
scheme uses three levels in the time direction. 



The conditions for the auxiliary variables Wa used in the Lindblom et al drivers are found 
from Q to be Wt(0,r) = Wr(0,r) = 0. 



5. Numerical Approach 

Here we describe our strategy for the numerical solution of the GH system (with a scalar 
matter source) in spherical symmetry. 

5.1 The numerical grid and the algorithm 

We cover the t-r plane by a discrete lattice denoted by (t^,ri) = {n At,i Ar), where n 
and i are integers and At and Ar define the grid spacings in the temporal and spatial 
directions, respectively. We note that when we perform convergence studies, we keep 
the ratio At/Ar constant so that our numerical scheme is generally characterized by a 
single discretization scale, h, which we can conveniently identify with Ar. As described 
in the next section, the spatial domain is compactified, and hence a grid of finite size A^^ 
extends from the origin to spatial infinity. As depicted in Fig. |^, approximations to the 
dynamical fields, collectively denoted here by Y, are evaluated at each grid point, yielding 
the discrete unknowns 1^" = y(t"',rj) = y(nAt,i Ar). In the interior of the domain, the 
GH equations and the gauge-driver equations are almost always discretized using 0{h'^) 
finite difference approximations (FDAs), which replace continuous derivatives with the 
discrete counterparts given in (|C.1|) and (|C.2|) . As in |^, |9| our scheme directly integrates 
the second-order-in-time equations (i.e. we do not rewrite the equations as a system which 
is first order in time). 

Following discretization, we thus obtain finite difference equations at every mesh point 
for each dynamical variable. Denoting any single such equation as 

Cy\7 = 0- (5.1) 
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we then iteratively solve the entire system of algebraic equations as follows. 

First, we note that for those variables that are governed by equations of motion that 
are second order in time, our 0{h?') discretization of the equations of motion results in 
a three level scheme which couples advanced-time unknowns at t"'^^ to known values at 
retarded times and In order to determine the advanced-time values for such 

variables, we employ a point-wise Newton-Gauss-Seidel scheme: starting with a guess for 
y^"^^ (typically, we take Yl^^^ = yj") we update the unknown using 



Here, TZy is the residual of the finite-difference equation ( |5.l| ), evaluated using the current 
approximation to Y^~^^ , and the diagonal Jacobian element is defined by 



= (5.3) 
dY 



In the cases where we used gauge drivers that involve B and Wa, we found that an it- 
eration based on an implicit Euler discretization scheme of the corresponding first or- 
der equations performed well.^'^ Specifically, writing any such equation schematically as 
Y = fyiY, dY, . . . ), we update using 

yn+l ^ yn-1 ^ ^At (5.4) 

We iterate { f).2\ ) and ( |5.4D over all equations until the overall residual norm^^ falls below 
some specified convergence threshold. 

In order to inhibit high-frequency^^ instabilities which often plague finite difference 



equations such as ours, we add explicit numerical dissipation of the Kreiss-Oliger type |31| 
to our scheme. Following at every grid point and for each dynamical variable we make 
the replacement 

Yi,^Yi- exo di (5.5) 
at both the t^~^ and time-levels before updating the t""*"^ unknowns. Here, di is defined 

by 

di = ^ {Yi^2 - 4 y,„i + 6 - 4 y,+i + y^+a) . (5.6) 
lb 

and eKO is a positive parameter satisfying < eKO < 1 that controls the amount of 
dissipation. An extension of the dissipation to the boundaries as well as to the black hole 
excision surface (see Sec. |5]^ ), was also tried, but was not found to have any positive effect. 
In fact, using dissipation at the outer boundary usually resulted in late-time instabilities 



^''The advantage of the implicit Euler method is that it is unconditionally stable and easy to implement. 
Although it is only first-order accurate — which does impact the overall convergence of the scheme when the 
Lindblom et al drivers are used — we have found it useful to achieve our chief current goal of constructing 
stable numerical implementations for our GH system. 

^^defined, e.g. as a sum of absolute values of the individual residuals of the equations, TZ — X]y l^y I- 
"High-frequency" refers to modes having a wavelength of order of the mesh spacing, h. 



-17- 



in the code. 



5.2 Coordinates and boundary conditions 

While the physical, asymptotically flat spacetime extends to spatial inflnity, in a numeri- 
cal code one can only use grids of finite size. A standard strategy to deal with this issue 
involves truncating the solution domain by introducing an outer boundary at some finite 
radius where approximate boundary conditions are imposed. When such an approach is 
adopted, it is then important to ensure that the computed solutions do not depend sensi- 
tively on the truncation radius. However, another technique which has been successfully 
used in previous work in numerical relativity, see e.g. ||9|, involves compactification of 
the spatial domain. Paralleling the experience of these earlier studies, we have found that 
compactifying the radial direction and imposing the (exact) Dirichlet conditions ( [4.1[1| ) at 
the edge of the domain works well, provided that we use sufficient dissipation. In partic- 
ular, it is known that due to the loss of resolution near the compactified outer boundary 
(assuming a fixed mesh spacing in the compactified coordinate), outgoing waves generated 
by the dynamics in the interior will be partially refiected as they propagate towards the 
edge of the computational domain, and these reflections will then to tend to corrupt the 
interior solution. By adding sufficient dissipation one can damp the waves in the outer 
region, attenuating any unphysical influx of radiation, and thus enabling a meaningful use 
of compactification. 

For the general case where we have more than one spatial dimension, X*, requiring 
compactification, we consider a transformation that maps X* G [0,cxd) onto x* G [0, 1], 

X' = O(x^), (5.7) 

where the Q are monotonic functions, such that Ci(0) = 1, and which will have essential 
singularities at = 1. The field equations ( [4.12 - |4rT^ are discretized in the compactified 



coordinates after we analytically remove the Jacobian of the transformation (qJ) in all the 
differential operators. The general replacement rule for first and second spatial derivatives 
is dx = eidx and dj^ = ef9^ -|- e2dx, where ei = !/(' and 62 = —C'/iCY^ so, for example, 
a typical term in ( 4.12 - 4T^) , dgu/dX^ , would be replaced with {(j)~^dgti/dxK 



In the spherically-symmetric calculations considered in this paper we use a specific 
compactification 

r = (5.8) 

I + r 

where the compactified f ranges from to 1 for values of the original radial coordinate 
r G [0, cx)). The boundary conditions at f = 1 are then imposed exactly: gu = —l,gtr = 
0,fi'rr = 1, A = S = 0, and = 0. For the gauge source functions we set Ha = 0, as well as 
Wa = B = 0. 

We have previously described the boundary (regularity) conditions at f = r = in 
Sec. L2. Denoting by Y^~^^ the advanced-time value at the origin for any of the variables, 
9tti9rr and Ht that have vanishing derivative at r = 0, we use the update Y^^^ = (4 1^2"^^ ~ 



Y^+')/2,, which is based on an 0{h ) backwards difference approximation (see ( p.3| )) of 
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drY = dfY = 0. For the quantities gtr and H^-^ which are odd in r as r — )■ 0, we simply use 



= 0. 



As discussed in Sec. [4.2| , we considered the introduction of a new variable, A ( 4.11| ), 



to expedite implementation of the regularity conditions involving and S. However, in 
the calculations described below we have adopted a simple method that does not involve 
A and that works well in spherical symmetry. In this approach, we retain the original 
variables S and grr, and impose = and S = (1/2) log(g'rr) at the origin. Then instead 
of determining (i.e. the advanced value of S at the next-to-extremal grid point) from 
the corresponding discrete evolution equation, we perform the update using the 0{h?') 
backwards FDA to the regularity condition, S'{t, 0) = 0, namely 5^+^ = (3 5^'+^ + 5J+^)/4. 

We must also maintain regularity at the origin for the auxiliary functions Wa and B 
that are used with some of the gauge driver conditions. We expand the metric functions 
in analytic Taylor series around r = and substitute the expansions into the equations 
( |0|J4.24| ) to arrive at 



^ + r?2 5 = 0, 
Wt + gtt [t] Wt grr -{n + l)H';) = 0, 

Wr + gtt {v Wr grr " H';) = 0, (5.9) 

which we use to advance B{t,0) and Wa{t,0) forward in time. Operationally, the time- 
derivatives in the equations are replaced with the FDA expressions ( |C.1D evaluated at t", 
and the spatial derivatives are replaced with one-sided versions ( |C.3| ) evaluated at 
The values of the functions B{t"'~^^,0) and Wa{t"^~^^,0) are then algebraically found. 



5.3 Apparent horizon and excision 

As is well known from many theoretical studies (both closed- form and numerical), a grav- 
itational collapse process that concentrates sufficient mass-energy within a small enough 
volume can lead to the formation of a black hole. In numerical calculations based on a 
space-plus-time split, black hole formation is often inferred by the appearance of appar- 
ent horizons. We recall that an apparent horizon is defined as the outermost marginally 
trapped surface, and that a marginally trapped surface is one on which future-directed null 
geodesies have zero divergence. Specifically, given a surface with outward-pointing space- 
like unit normal, s", embedded in a hypersurface with future-directed timelike unit normal, 
71°, the vanishing of the divergence, 9, of the outgoing null rays defined by P = + n" 
can be expressed as 

= _ s-/)V„/;3 = 0. (5.10) 



However, we have checked that the scheme that uses A performs remarkably well in our 2 + 1 numerical 
implementation [P2| that generalizes the present 1 + 1 work. 
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In spherical symmetry we have = grr^"^ dr, and the above equation can be written as^^ 
9 = rdtS + {l + rdrS) I + J ^ - ^] =0, (5.11) 



In numerical calculations, one can thus easily locate an apparent horizon by simply search- 
ing for zeros of 6: the position of the outermost such zero then coincides with the location, 
TAHi of the apparent horizon. 

In our code we use excision to (dynamically) exclude from the computational domain 
a region interior to the apparent horizon that would eventually contain the black hole sin- 
gularity. The success of this approach hinges on the observation that in spacetimes that 
satisfy the null energy condition (such as those that we construct) and assuming cosmic 
censorship, the apparent horizon is contained within the event horizon, which ensures that 
the excluded region is causally disconnected from the non-excised portion of the domain 
(see [^] and the references therein for further discussion). Operationally, once an ap- 
parent horizon is found, we introduce an excision radius, texi that satisfies tex < ''AHi 
and such that all radial characteristics at r = tex are pointing inwards. (We typically 
find Tex ~ 0.4 taH) where we again emphasize that r is the coordinate radius.) This spe- 
cific characteristic structure eliminates the need for boundary conditions at tex^ rather, 
advanced-time unknowns located on the excision surface are computed using finite dif- 
ference approximations to the interior evolution equations, but where centered difference 
formulae are replaced with the appropriate one-sided expressions given by ([ 



5.4 Spacetime diagnostics 

We employ several diagnostics in order to characterize the geometries of the spacetimes we 
construct. 

Mass. Far away from an isolated system a natural radial coordinate is defined by the 
asymptotic flatness of the spacetime, and the ADM mass of the solution can be found from 
the asymptotic radial behavior of the metric functions. In spherical symmetry there is only 
one asymptotic constant, tq, that can be determined, for instance, from the fall-off of gu- 
gtt ^ ^ + rQ^^/r''-^. This constant is related to the mass ^ hy M = nfin/(167r)ro"\ 
where iln = 27r*^^""*'^^/^^/r[(n + l)/2] is the surface area of a unit n-sphere. 

In addition, in spherical symmetry one can define a local mass function, m{t, r), some- 
times called the mass aspect 

m{r, t) = ^^^j^ (l - RaRpg''^) , (5.12) 

where R = re^ is the areal radius. The mass aspect is negative inside a trapped (or 
anti-trapped) region, vanishes at its boundaries and is positive outside in regular region. 



^*An alternative way to derive this result relies on the fact that the apparent horizon in spherical symmetry 
can be defined as a null surface located at constant radius. Equating the time-derivative of the areal radius 
along null rays to zero, d{r e'"^) /dt\ia, = re^dtS + r (1/r + drS) ^-gtr/Srr + \/ff?r/5rr - 9tt/grr^ = 0, 



where the expression in the second brackets is dr/dt\ic, we recover the result in (5.11) 
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Figure 2: Outgoing null rays in the t — R plane emitted from the origin at different times {R is 
the compactified areal radius). The left panel shows the geometry generated by an initially origin- 
centered pulse of matter with $o = 1-6 that disperses infinity. The presence of the matter deflects 
the outgoing null rays towards the origin, but the rays eventually escape to infinity. The motion of 
the pulse can clearly be traced. The right panel shows the geometry generated by stronger initial 
data having $o = 3.0. In this case the matter collapses to form a black hole of mass, A/bh — 0.3: 
rays emitted before t ~ 1.25Mbh escape to infinity but the rays emitted after that time fall back to 
the origin. The null ray that separates the two regimes designates the event horizon and the thick 
dashed line is the asymptotic apparent horizon. The thin dashed lines are obtained by integrating 
(5.13) backward in time, and are attracted to the event horizon. 



It grows monotonically and asymptotically coincides with the ADM mass. 

Null geodesies. A convenient way to visualize the causal structure of a spherically 
symmetric spacetime is to plot a family of outgoing null rays, la- When plotted in the t-R 
plane, the slope, dR/dt\ia, of an outgoing null geodesic is positive outside the apparent 
horizon, and asymptotes to the flat-space value of unity for large values of R. Additionally, 
the slope vanishes at the apparent horizon, concomitant with the vanishing of the outgoing 
null divergence, and becomes negative inside the horizon. All of these features can be seen 
in Fig. where the displayed lines are integral curves, i?(t;to)- Here R is the compactified 
areal radius, and the corresponding uncompactified trajectory, R{t;to), is defined by 



R{t-to) 



to 




OR dR 
dr dt 



dt'. 



(5.13) 



Each curve thus represents the path of an outgoing null ray that is emitted from the origin 
at a specific time, i = to- 
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Event horizon. In contrast to the local definition ( 5.11 ) of the apparent horizon, the 
event horizon is a global concept: it is defined by outgoing null rays that neither escape to 
future null infinity, nor fall into the black-hole singularity. Clearly, this definition requires 
knowledge of the complete time evolution of the system, and hence, assuming a calculation 
that is carried out for a finite amount of coordinate (or proper) time, one cannot even 
in principle locate event horizons in numerically-generated spacetimes. However, when 
a spacetime approaches a stationary state, an approximate event horizon can be found. 
We employ the method of Libson et al (sj] which is based on the observation that if one 
integrates the geodesic equation ( 5.13| ) backward in time, the event horizon becomes an 



attractor for geodesies that either escape to future null infinity or fall into the singularity 
at arbitrarily late times. We have found that in our simulations the event horizon is traced 
fairly well by the time development of the apparent horizon. Again this can be seen in 
Fig. H, where the thin dashed lines show the trajectories obtained by integrating ( 5.13| ) 
backwards in time, and starting with several initial radii. 

6. Results 

For concreteness, we restrict our numerical experiments to the case of four-dimensional 
spacetimes, and take our matter source to be a real, massless scalar field. All of the 
results discussed here were generated using an initial scalar field profile of the Gaussian 
form ( 4.25| ), with fixed values tq = and A = 0.6, so that the scalar pulse is always 



initially centered at the origin. The overall amplitude, of the profile ( f4.25| ) is then used 
as a control parameter: variations of produce varying "strengths" of initial data, and 
varying degrees of non-linearity in the ensuing evolution. In practice, the maximum value 
of 2m{t,r)/R(t,r) (where R is the uncompactified areal radius) that is achieved in a given 
calculation is a useful indication of how strong-field the evolution becomes. 

We use the above notion of initial data strength to loosely define three classes of 
solutions — within a given class we observe that the overall dynamics of each of the scalar 
and gravitational fields are similar. Specifically, we consider the following cases: (i) weak 
data, defined by $o ^ 0.5, yielding maxi_r2m/i? ~ 0.08; (ii) intermediate data, having 
0-5 < ^0 ^ 1-6, and maxf^r2m/i2 ~ 0.25, and (iii) strong data, with $0 > 1.6 and 
maxt^r 2m/ R > 0.25. While the first two cases describe weakly and mildly gravitating scalar 
pulses, respectively, which completely disperse in all instances, the strong data generates 
spacetimes in which black holes form, or almost form (i.e. near-critical evolution, see (p7|])). 

We have also found it useful to use the total ADM mass, Madm, of the spacetime — 
which can be computed at t = — to normalize certain numerical parameters. In particular, 
we set the parameters of the gauge driver ( 3^ , 3^ ) using ^1 = Cio/-^adM' ^2 = '^2o/-^adm 



and ^3 = '?3o/-^ADM' where the "bare" values, kq, ^10, C20 and ^30 are generally held fixed 
as ^0 is varied. Moreover, and as discussed in more detail below, we find that the accuracy 
of our results is improved if the constraint damping term asymptotically vanishes at large 
spatial distances. Accordingly, we typically multiplied k by the factor 2Madm/-R- 

Because we use, at least in large part, a time-explicit finite difference scheme, we 
expect restrictions on the ratio Ac = At/Ar (the Courant factor) that can be used while 
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maintaining numerical stability. For the case of harmonic gauge, we found that values of 
Ac satisfying 0.01 < Ac < 0.8 generated stable solutions with roughly constant accuracy, 
although somewhat stronger numerical dissipation was required to stabilize runs that used 
larger values of Ac in that interval. In the results discussed below we have typically taken 
0-3 ^ Ac ^ 0.6 for weak and intermediate data, and 0.1 ^ Ac ^ 0.2 for the evolution of 
strong data. We further found that when any of the other gauge drivers were adopted, 
smaller Courant factors (relative to the harmonic case) were required. In those cases our 
results were generally computed using 0.05 ^ Ac ^ 0.1. Typically, in cases where Ac 
was taken too large, we observed amplification and dominance of numerical errors near 
the origin: this lead to high frequency oscillations and, eventually, to divergence of the 
numerical solution. 

Another crucial numerical parameter is the Kreiss-Oliver dissipation factor, exO) which 
we generally set according to 0.1 < exo ^ 0.7. Finally, it is important to note that we 
found that optimal values of both Ac and exo were dependent on the spatial resolution: 
specifically, as Ar — )• somewhat smaller values of Ac, as well as larger values of exo were 
usually required. The lowest and highest resolution runs reported in this paper typically 
had Ar =1/64 and Ar = 1/8192, respectively: runs with Ar = 1/8192 generally required 
Ac = 0.05 and ^ko = 0.7 for stability. 

Many of the coordinate conditions discussed and employed in this paper are charac- 
terized by several adjustable parameters, and we have by no means carried out exhaustive 
parameter space surveys in all cases in an attempt to optimize parameter settings. Rather, 
our more limited numerical experimentation indicates that with a certain amount of tuning 
of the parameters, it does seem possible, at least in principle, to simulate various interest- 
ing situations. Our intent here is chiefly to document the overall behavior of several gauge 
conditions as well as to explore some of the effects that specific parameters of the gauge 
drivers have on the evolution. Given this primary goal, we also defer most of our discussion 



of code convergence and accuracy to Sec. pA . 
6.1 Weak data 

In this section we consider the evolution of weak initial data for which $o ^ 0.5, yielding 
-^ADM ^ 0.01 and maxj^^ 2m/ R < 0.08. In this case there is little interaction between the 
scalar and gravitational fields, the scalar pulse entirely disperses to infinity, and we find that 
essentially any of the gauge conditions described above can be used to produce long-term 
stable evolution. For this type of data we use exo — 0.1 for the Kreiss-Oliger dissipation 
parameter, finding that larger values have detrimental consequences for stability. However, 
even with dissipation and constraint damping, we find that numerical errors eventually do 
grow — on a time scale of order t > IO^Madm — and cause the code to crash. 

We find that the effect of the constraint damping term depends on whether k is fixed 
or allowed to vary over the integration domain. For fixed k, it is essential to take kq > 0.01, 
otherwise high-frequency oscillations quickly ruin convergence. However, if the damping 
is too strong, instabilities are also triggered. In fact, we find that the optimal damping 
parameter is related to the typical scale over which the scalar field varies. For the Gaussian 
initial data that we consider, this scale is A, so we take k ~ A~^. (This observation holds 
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for intermediate strength data as well, as can be seen in Fig. ^.) On the other hand, when 
we take k = K{r), and specifically for the choice k = ko{2Mabm/ R) mentioned previously, 
we find that the results are relatively insensitive to the value of kq, provided kq ^ 100A~^. 
For larger values of kq instability is again usually observed. 

Our experiments with the gauge drivers proposed by Lindblom et al, have focused on 
the specific Bona-Masso slicing condition for which /(a) = 2/a, corresponding to 1 + log 
slicing. However, for weak data, we find that other choices of / (such as f{a) = 2a, o? and 
10/q^, to list a few that we have tried) produce qualitatively similar results. 

Considering the conditions that determine the shift, we find that the F-driver condition 
performs somewhat better than F-freezing, with the former allowing the evolution to be 
controlled for a longer amount of time. There was only mild dependence on the gauge- 
driver parameters, fii^2,Vi,2,Qs, gn,(^ and provided they are all taken in the range 0.01-10 
in units of Madm- 

In order to assess the performance of the coordinate conditions in driving the source 
functions to the target functions, we first follow [^] and define the weighted L2-norm, \Y\, 
of a function Y as follows^^. 



A similar, if somewhat less smooth norm, which we also use here, can be defined as 



Nr 



i=l 



Fig. ^ shows the weighted norms of the differences between the actual and target 
source functions from a typical weak- field simulation. It is evident from these plots that 
the drivers successfully drive the source functions Ha towards the target functions Fa as 
the evolution proceeds. 

We now continue to discussions of the evolution of intermediate- and strong-field data, 
where the results are more sensitive to the specific driver used, as well as to the parameter 
settings for any given driver. 

6.2 Intermediate data 

Here we consider evolutions characterized by 0.5 < <I>o ^ 1-6, where Madm ^0.1 and 
maxt^r 2m/ R < 0.25. First, for this strength of data, we have found that the pure harmonic 
and GH gauges ( |4.15 - 4T6| ) perform comparably. With both choices, we are typically able 



to accurately trace the evolution of the initial data for times of the order of 100-600 Madm, 
with increasing resolution resulting in increased maximum evolution time. 

The causal structure of the spacetime from a typical intermediate strength computation 
is displayed in the left panel of Fig. We recall that in this figure the curves represent 



^The integrals are evaluated on our fixed mesh using the trapezoidal rule. 
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Figure 3: The behavior of gauge drivers for the case of Bona-Masso shcing with /(a) = 21 a (left) 
and the F-driver shift condition (right) in the weak field regime, $o — 0.1 



trajectories of outgoing null rays that are emitted at regular intervals (in coordinate time) 
from r = 0. As the evolution proceeds, the pulse, which is initially centered at the origin, 
disperses to infinity. The outgoing null rays are bent towards the origin by the presence 
of the matter and asymptotically become straight lines with unit slope in the r — t plane. 
The position of the scattered pulse of scalar field can be traced through the location of 
the "ripple" in each curve, i.e. at the positions where the outgoing null geodesies suffer the 
most deflection. 

We will discuss issues of code convergence and accuracy in more detail in Sec. |6.4| . 
However, we note here that constraint norms, IMqI^j) defined by ( B.ll| ) and computed, 
for example, using either ( |6.1| ) or ( |6.2| ) provide a basic indication of the accuracy of our 
numerical method. For the calculation depicted in Fig. |2| that uses a medium resolution, 
Ar = 1/1024, we find the initial norms iM^li^jOf order 10~^, which for roughly the first 
half of the evolution then decrease to values of 10~^--10~^. Thereafter we observe a slow 
increase in the size of the constraints although — except for the last few time steps before 
the code fails — |Mq,|2,2 remain well below the 10~^ level. Moreover, we generally observe 
the expected quadratic convergence of iM^lij as the finite difference mesh is refined. 

Another basic indication of numerical accuracy is provided by the the sum of the norms 
defined by ( |6.2| ) of the residuals of the dynamical equations, \R\l\ = X^y l^^'l' where TZy 
is the FDA residual of the equation that governs the field Y . Fig. ^ shows the behavior 
of |7^| function of the damping parameter, kqi for calculations with = 513 

(moderate resolution), <I>o = 1.6, and where k = ko(2Madm/^)- As already noted in the 
discussion of the weak field results, the sizes of the constraint and equation residuals tend 
to be minimized when kq is comparable to the inverse of the typical length scale of the 
problem, i.e. to for our initially Gaussian data. This is apparent in the figure, which 
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Figure 4: Illustration that the characteristic behavior of the L2 residuals of the evolution equa- 
tions depends on the value of the damping parameter. Excessive or insufficient damping degrades 
convergence, or leads to divergence. The optimal range for the damping parameter is kq 
where A is the typical length scale in the problem. 

shows that for kq = 0.5/A, the residuals remain on the order of 10~^. 

We next experiment with the Lindblom et al drivers, and find that while for $0 < 0.7 
the dynamics of Ha and is qualitatively similar to that in the weak field regime (shown 
in Fig. ^) and essentially independent of the parameters of the gauge drivers, for $0 > 0.7 
the convergence of the source functions. Ha, towards the target sources, Fa, has stronger 
dependence on the parameter settings. The most pronounced feature in this regime is that 
the drivers succeed in forcing Ha — )• Fa only on the length-scale set by the parameter //i. 
In particular, when we start with initial data that has Ha = Fa, we find that for large 
values of ^1 the source functions remain close to their targets for a a few tens of MadM; 
after which high-frequency oscillations destroy the matching. Conversely, starting from the 
same initial set up, but taking ^1 very small, we observe that the source functions quickly 
deviate from the targets and never approach them in the subsequent evolution. 

Given this observation, and given that our Gaussian initial data generates an evo- 
lution characterized by a length scale, A, it is thus reasonable to take /ii ~ 1/A in an 
attempt to enforce the desired gauge conditions on that scale. Results from such a com- 
putation are shown in Fig. ^, which displays the source and target functions, as well as 
their Fourier transforms, from the evolution of initial data with ^>o = 0.9. The calculations 
were performed using target slicing of the Bona-Masso type with f{a) = 2/a, and target 
F-driver shift conditions with /^i = 1.3 (recall that A = 0.6 for all of the computations de- 
scribed here). In addition, here, and for all of the results discussed in this section, we used 
IJ^2 = V = ^2 = ^, Qn = Qs = 0.5,(7 = —1/3 and = 0.7. In contrast to the case of ni, we 
find that the calculations are not too sensitive to the settings of these parameters, so long 
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Figure 5: The source functions i?", the target functions F", and their Fourier transforms at two 
instants. We begin with initial data satisfying Ha — Fa. Within a few dynamical times the functions 
deviate, but subsequently are driven towards each other. After a time of 30 — 50A/adm they match 
on length-scales of order This is illustrated by the spatial spectral decomposition shown in 

the right panels: while the lower frequencies of the functions match closely, the higher-frequency 
components do not. 



as their values are all of order unity. In this simulation we begin with initial data satisfying 
Ha = Fa. Within a few dynamical times the functions deviate, but as Fig. || demonstrates 
the functions are subsequently driven towards each other, when the source functions start 
resembling the targets on the spatial scales Notice that the high-frequency spatial 

variations of the target F°'s are not replicated by the source functions. Similar behavior 
was originally observed in |18] for perturbations on a given background. 

The manner in which the coordinate conditions evolve in time for this calculation is 
shown in Fig. which depicts the norms of the functions Gt and Gr, defined by ( |3.9| ) 
or ( 3.1C| ), and ( |3.14| ) or ( 3.17 ). As described in Sec. |3|, enforcing a particular gauge is 
equivalent to driving these functions to zero. Since we begin with initial conditions in 
which the gauge is exactly fixed, the norms of Gt and Gr are initially zero. Then on a 
timescale of order several tens of MadM; the norms grow to some maximum value, after 
which they decrease slowly. The details depend on the particular coordinate choices, as 
well as on the settings of the driver parameters, but usually it is possible to drive the 
L2-norms of Gt and Gr to the level of about 0.01. 

Although for smaller initial pulse amplitudes ($0 ^ 1-0) we managed to find parameters 
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Figure 6: The norms of the functions Gt and Gr, that must vanish when a specific gauge to which 
they correspond is approached. The norms decrease in time in a way that depends on the details of 
the gauges and the parameters of the drivers. In the simulations shown, we mostly use Bona-Masso 
slicing with various choices of /(a), and F-driver shift conditions, except for the data plotted with 
+-symbols that was obtained using a F-freezing shift condition. In all cases we are able to drive 
the norms to a level of about 0.01. 



for the Lindblom et al drivers that asymptotically fix the desired gauges, we find that for 
larger amplitudes the effectiveness of the drivers degrades, and for <I>o ^ 1.0 we could not 
find parameter settings that enforce any of the specific gauges. This does not necessarily 
mean that the code diverges: indeed, the evolution often proceeds, but the behavior of the 
source function is rather arbitrary. In this regime we find that the evolution systems based 
on the Lindblom et al drivers tend to be more dynamical and less stable than one that uses 
simple drivers such as (4.15). 



6.3 Strong data and black hole formation 

Increasing the initial amplitude, ^q, of the scalar pulse leads to increasingly strong cur- 
vature in the development of the initial data. As expected, above a critical value — in the 
current case, <l>o ~ 2.15 — black holes form, as signaled by the appearance of apparent 
horizons. We recall that we have already used the trajectories of outgoing null geodesies 
to schematically display the causal structure of a typical black hole geometry in the right 
panel of Fig. ^ 

Our first set of numerical experiments in the strong-field regime compares subcritical 
evolution (<l>o ^ 2.15) in pure harmonic coordinates to that in the generalized harmonic 
gauge given by ( 4.15| ). A generic feature of purely harmonic evolution in this case is a fairly 



quick collapse of the lapse function towards zero values near and at r = 0. As a result 
the evolution in the central region (where the pulse is concentrated) effectively freezes, 
and the scalar field remains present near r = even at late (coordinate) times. This is 
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Figure 7: Proper tirae (6.3) at the origin in harmonic evolution as a function of coordinate time for 
several initial data strengths. The evolution slows down for stronger data and it effectively freezes 
for near critical data. 



demonstrated in Fig. ^, which shows the evolution of central proper time 

r(t) = [ a{t',0)dt', (6.3) 







function of the strength of the initial data. 

On the other hand, and in accordance with the previous experience of Pretorius |9|, 
we are able to use the generalized harmonic gauge condition ( [4.15| ) to inhibit the collapsing 
of the lapse. Specifically, we use ao = 1 and g = 3 in ( [4. IS]) , and experiment with 
various values for and ^2- In addition, motivated by an observation that we can more 
stably evolve subcritical data by gradually "turning-off" the gauge driving at late times, we 
actually replace and 6 in ( [ilsD by (eio/M2^^)/(l + stP) and (6o/Madm)/(1 + stP), 



respectively, where p and s are additional positive parameters. In practice, we have usually 
taken p = 1, leaving s free to control the rate at which the gauge driving is disengaged. 

Results from calculations with = 1.8 (Madm — 0.125) and using several sets of 
values for ^lo, ^20 and s are shown in Fig. |^. The plots clearly show how judicious choice 
of the parameters can prevent the collapse of the lapse. Through experiments with various 
subcritical initial data sets we find that parameter values 1 < ^10 < 5 and 0.5 < ^20 ^ 2 
produce good results. However, in order to keep the lapse from collapsing for initial data 
very close to criticality, we generally needed to increase both ^10 and ^20 by factors of 
up to 10, while simultaneously increasing s (to values of order 50) and taking p = 2 or 
3. For instance, simulations that use 2049 spatial grid points and the driver ( 4.15| ) with 



the parameters tuned to ^10 = 50,^20 = 30, s = 36 and p = 2 allowed us to explore the 
dynamics of solutions with = 2.1465 it 0.0005 without encountering a collapsing lapse. 
Unfortunately this is not close enough to the threshold amplitude for us to be able to 
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Figure 8: Left panel: the proper time (6.3) at r = from an evolution tliat uses gauge conditions 
(4.15) with $0 = 1-8, Hr = 0, g — 3, and and ^2 additionally divided by 1 + s i. Right panel: the 
amplitude of the scalar field at the origin from the same simulations. While in harmonic gauge the 
evolution freezes near r = 0, in the dynamical gauge (4.15) it continues. 



observe in detail the distinctive features of scaling and echoing known to appear in the 
near-critical regime of this model |37]. 

We end our discussion of subcritical strong-field evolution with two observations. First, 
we note that while we have investigated the use of dynamical conditions such as ( 4.1(^ ) 
for Hr, the spatially harmonic choice, Hr = 0, is simpler to implement, and apparently 
more stable in this regime. Secondly, although we have experimented extensively with the 
Lindblom et al drivers in this context, we have not been able to find parameter settings that 
prevent coordinate pathologies (premature collapse of the lapse) from quickly developing 
for near-critical evolutions. 

We now turn to the case of supercritical evolutions, which are characterized by the 
formation of black holes. As described in Sec. we have implemented black hole excision 
techniques in our code: however, due to the strong singularity avoidance property of pure 
harmonic gauge, as well as the generalized harmonic modifications (4.15- 0!^ ), we can also 
perform computations in which black holes form and are evolved for some amount of time, 
but where excision is not used. 

For example. Fig. |^ shows metric functions from a calculation with <I>o = 3.0 that 
uses pure harmonic gauge with no excision. We infer the formation of a black hole by 
the appearance of an apparent horizon, which at the end of the simulation is located at 
a compactified areal radius i?AH — 0.4. We can then estimate the mass of the black hole 
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Figure 9: Illustration of the geometry of black hole formation without excision. The metric 
functions remain regular all the way to the origin, where the functions tend to zero. The black hole 
which forms has a mass Mbh — 0.335 with a horizon at i? ~ 0.4 in the compactified areal radial 
coordinate. Note that the lapse collapses at the origin, freezing the evolution there. 




5 10 15 20 25 



bh 

Figure 10: Time plot of the central value of the Kretschmann scalar, showing indefinite growth 
which signals the development of a curvature singularity. 

at that time from the apparent horizon location: Mbh = 0.5 -Rah — 0.34, and note that 
the total ADM mass in this case is Madm — 0.41. An apparent horizon is first detected 
at t ~ 1.25Mbh and Fig. ^ displays the metric functions at two instants: (i) t ~ 5Mbh 
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Figure 11: A strong data simulation, $o = 3. Shown are Z/2-iiorms ( pj| ) of the normalized 
deviations between the target and source functions (top), and the gauge functions, G", which 
vanish when the desired gauge is achieved (bottom). The source functions, Ha, approach the 
targets. Fa, uniformly soon after the horizon forms at t ~ 3 A/bh and follow them closely until low 
frequency variations in F", induced by the matter outside the black hole, develop and destroy the 
uniform match. The sources replicate the targets (not shown) on the scale defined by which 
we take to be of order of i?AH- 



(dashed lines), and (ii) t ~ 20MbH) which is shortly before the simulation crashes (solid 
lines). For this specific calculation we used 4097 spatial grid points, and, at the time of 
the code crash, the values of the temporal component of the metric, gu, near the origin 
are of order 10~^^ (corresponding to lapse values of order ~ 10~^). Despite the fact that 
all of the metric components displayed in Fig. ^ are tending towards zero at the origin 
at late times, the functions remain smooth and regular throughout the evolution. Fig. 10 
plots central values for the Kretschmann scalar, Ra^-ysR"'^'^^ , as a function of time. The 
apparent divergence of this geometric quantity indicates the development of a curvature 
singularity. 

We have found that the use of excision can somewhat extend the duration of our 
simulations of black hole spacetimes. For comparison, a run with the same parameters 
enumerated above, but employing excision, lasted for as long as ~ GOMbh- We recall that 
our simple approach to excision has been described in Sec. 5^, and note that in practice 
we have typically chosen the excision radius, tex, to satisfy tex ^ 0.4rAH- The rest of the 
results described in this section were obtained in simulations with excision. 

Although we are able to avoid the central singularity using excision, it is clear from 
our calculations that the harmonic coordinate system continues to evolve in a highly non- 
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trivial manner after excision is initiated. This dynamics in the coordinates causes, or is at 
least associated with, two main problems. First, the resulting coordinate system does not 
approach a stationary state: in particular, the coordinate position of the apparent horizon 
evolves with time. Specifically, after formation, the horizon expands outwards and con- 
sumes most of the numerical grid. Eventually then, the portion of the spacetime outside 
the horizon — which we recall extends to spatial infinity due to our use of a compactified 
coordinate system — is represented by only a small portion of the initial lattice. Conse- 
quently, numerical errors that arise near the outer boundary dominate the late stages of 
the evolution. The second problem is that the lapse continues to decrease in the vicinity 
of ^EXi and becomes very small. In this situation truncation errors in quantities near tex 
occasionally cause the computation of non-positive values for the lapse, which immediately 
leads to code failure. Both of these problems can be somewhat mitigated by increasing 
the numerical resolution. In harmonic gauge, we were able to simulate the formation of 
a black hole and resolve it for about t ~ 70Mbh using our finest resolution, A''^ = 8193. 
However, given these difficulties induced by the late-time dynamics when using harmonic 
coordinates, it is quite natural to try to use the coordinate freedom provided by the various 
gauge drivers discussed above to a) attempt to minimize the time development of the lapse 
following the formation of an apparent horizon, and/or b) implement a non-trivial shift 
vector with an aim to minimize the outward expansion of tah at late times when there is 
very little matter falling into the black hole. We thus now summarize our experimentation 
with several driver conditions that was focused on realizing these ideas. 

As we have already mentioned, one of the main motivations for Pretorius' development 
of the driver condition ( 4.15| ) was to keep the lapse from collapsing in the vicinity of 
horizons [^. Following that work then, we first used (|4.15| ) to fix the time slicing, while 
maintaining harmonic spatial coordinates {Hr = 0). However, in contrast to the results 
reported in (which we note were performed in three spatial dimensions using Cartesian 
coordinates), we found the evolution in this case to be significantly less stable than purely 
harmonic evolution. For example, even a small value of of order O.OlMgjj resulted in a 
code crash at a time about a factor of two earlier than for the harmonic case, irrespective 
of the value of the friction parameter, ^2- 

We next used harmonic slicing, Ht = 0, while evolving Hr using the driver ( 4.161 ). 
Here, we found a modest amount of improvement over the purely harmonic case, in that the 
"grid-sucking" phenomenon described above was slowed, with an accompanying reduction 
in the development of numerical error in the outer, low-resolution region. For example, the 
duration of the evolution of <I>o = 3 initial data with ^20)^30 ~ O{10)Mbh x (1 + 5t^)~^, 
increases by approximately 20% compared to the corresponding harmonic evolution. 

Interestingly, we obtained even better results using certain versions of the Lindblom et 
al gauge drivers. For the strong-field, supercritical calculations described here, we found 
that versions of the drivers that use the simple scalar operator ( [1.19 ) performed better than 
those that used ( 4. IS ). Moreover, we found that drivers based on the Bona-Masso slicing 
and F-driver shift conditions (with suitably tuned parameters) gave the best results, and 
for convenience will hereafter refer to this specific choice as BMGD. In particular, relative 
to other driver choices, this combination minimized — but unfortunately did not completely 
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Figure 12: Illustration that the level of constraint preservation is dependent on the choice of kq, 
and that the optimal value is in the range kq ~ . Excessive damping leads to rapid and violent 
growth of the constraints and divergence of the solution. 



eliminate — the outward drift of tah- Our best configuration allowed for accurate simula- 
tion of black hole spacetimes for about IOOMbh following the formation of an apparent 
horizon. After that time, code accuracy typically degraded, numerical errors near the ex- 
cision became dominant, and a late-time instability ensued. Based on our experiments, 
it remains unclear whether specific parameter choices for the drivers exist that would to- 
tally eliminate the drift of the coordinate position of the apparent horizon and, even more 
importantly, the disastrous collapse of the lapse inside the horizon. 

We now proceed to some details concerning our experience with the BMGD version 
of the Lindblom et al coordinate conditions. The parameters qn and qs that appear in the 
driver definitions — see equations ( 3.11 , 3.18 ) — control the relative weight that the gauge 
functions, G", have in forming the target sources, . We also recall that the G° vanish 
when the specific gauge to which they correspond is attained. We found it crucial not 
to choose qn too large: usually values in the range 0.01 — 0.1 resulted in the most stable 
evolutions, and would eventually lead to the desired behavior, — )• and G — t- 0. Our 
implementation was less sensitive to the value of qs, with results of comparable accuracy 
and stability being attained for Qs in the range 0.01 — 10. 

Having determined good values for qn and qs, we found through further experimenta- 
tion that stability is improved when the parameters ^1,^2 and rj are multiplied by a decay 
factor 2Mat)m/R in the region external to the horizon. This localizes the effect of the 
coordinate drivers to the near-horizon region, while producing a smooth blend to harmonic 
coordinates at spatial infinity. In addition, and in analogy to what we did for the subcrit- 
ical calculations in generalized harmonic coordinates described earlier in this section, we 
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further scale fii, 1x2 and 77, as well as g„ by 1/(1 + sf^). Here, s and p are again positive 
tunable quantities — we typically used p = 2 and s = 5 — that result in a late-time decay 
of the scaled driver parameters. We note that the values quoted below generally refer to 
"bare" values for parameters, with the additional scaling factors being implied. 



Fig. 11 shows the time development of the deviation between the target and actual 
source functions, Fa and Ha, respectively, as well as the gauge functions, G", for a typical 
BMGD calculation. The computation was performed with /xi = 4 ~ 1/Mbh, 1^2 = "T] = 
r/2 = 10 = 1, 0" = — 1/3,(?„ = 0.1 and qs = 1. The behavior of the two upper plots in 
the figure reflect the fact that the H'^ tend to the target source functions soon after an 
apparent horizon forms. Detailed examination of the data reveals that the match between 
the target and actual source functions is good throughout the entire domain for a certain 
amount of time following horizon formation. At late times the level of global agreement 
degrades, due to large scale variations in the F°' induced by the portion of the scalar field 
that is scattered to infinity. Despite this, we still find that actual sources accurately match 
the targets on the scale defined by /x^^^ ~ -Rah (not shown). The plots of the L2 norms of 
the gauge functions, G", shown in the bottom half of the figure, reveal a steady decrease 
in time, signaling that the desired gauge is being approached asymptotically. 

Our investigations of versions of the drivers using target functions corresponding to 
"static" gauges, such as maximal slicing and F-freezing, were unsuccessful in the sense 
that we were not able to find parameter settings that resulted in G" — )• as i — )• cx). 
Interestingly, however, we found that black holes could nonetheless be simulated using these 
conditions, with observed stability properties similar to those obtained using "dynamic" 
gauge conditions such as BMGD. This indicates that, at least for the type of initial data 



considered here, the stability of the drivers (3^,3^) does not strongly depend on the target 
gauge. 

Finally we note that the use of an appropriate amount of constraint damping is im- 
portant for computations in which black holes form. Fig. |l^ shows the behavior of the 
sum of the L2-norms of the constraints, ( |6.2D , in a sample run with N^. = 4097 and using 
various values for the damping parameter, kq. The plots provide clear evidence that the 
level of constraint maintenance (as well as the maximum simulation time) is optimized 
for Ko — -^BH- Values of kq significantly larger than the optimal value produce rapid 
code crashes, while those that are significantly smaller lead to poorer preservation of the 
constraints. 

6.4 Code accuracy, convergence and constraints 

In this section we briefly discuss some of the technical issues relating to the basic perfor- 
mance of our numerical code, including resolution requirements and checks of convergence. 

Not surprisingly, we find that the minimum discretization scale required to produce 
an acceptable evolution (for fixed choice of coordinate conditions) depends on the strength 
of the initial data. For example, in the case of weak and intermediate initial data, as 
defined previously, even a modest lattice size of A',. = 65 is enough to allow for long-time 
evolution. However, for stronger data, meshes sizes of at least N^- = 257 are required. 
Additionally, our code cannot evolve strong-field data for arbitrary amounts of coordinate 
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t/M„^, =120 

ADM 




Figure 13: Plot showing tliat tiie convergence of tlie scalar field is second order over most of the 
domain, with some irregularities occurring near the outer boundary and at the location of the scalar 
pulse. In this simulation, $o = 0.55. 




Figure 14: A log-log plot of the _L2-norms of the Hamiltonian constraint, Mt, and the momentum 
constraint, Mr for 5 different grid resolutions, and with $o = 0.55. The constraints remain small 
during most of the time-evolution, except for the last moments of the simulation, when instabilities 
set in and eventually lead to code failure. For the most part, the constraints converge as the 
resolution is increased, but there is a slowing of convergence at the highest resolutions. This issue 
is still unresolved, but may be related to the nature of the time iteration. 



time: generically, numerical problems develop that lead to a code crash on the order of 
10-100 Madm, and the precise lifetime of the simulation is dependent on the strength of 
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Figure 15: Linear plot of the L2-norm of the Hamiltonian constraint for 5 different resolutions 
(the same as those used in Fig. |4|), and again with $0 = 0.55. Modulo the comment in the caption 
of Fig. second order convergence is observed. 



the initial data, the resolution, and the details of the coordinate conditions. 

Much of the build-up of error that eventually leads to code failure, especially in sub- 
critical simulations, can be traced to the use of spatial compactification. In all of our 
calculations, there is outflux of scalar field to spatial infinity, and as the scalar radiation 
propagates to large distances it becomes more poorly resolved on the mesh, which has uni- 
form spacing in the compactified radial coordinate. Untreated, this will lead to spurious 
reflection of the waves which will corrupt the interior solution, so we add Kreiss-Oliger dis- 
sipation to explicitly damp the radiation when its wavelength becomes of order the mesh 
scale. Although this damping is imperfect, we find that increasing the resolution is effective 
in extending the lifetime of our evolutions. As a specific example, for a calculation which 
forms a black hole of size -Rbh — 0.6, and that uses BMGD coordinate conditions and 
excision, a grid with A^^ = 4097 is sufficient to keep the reflections small during all stages 
of the evolution until t ~ IOOMbh- Thereafter, an instability appears near tex and leads 
to a code crash. 

A crucial test of any finite difference code for the solution of a system of partial differ- 
ential equations involves the investigation of the convergence of the generated numerical 
solutions as a function of resolution. We perform straightforward convergence tests based 



on the assumption (originally due to Richardson |39|) that for any of the unknown functions, 
Y{t,r), appearing in our differential system, the corresponding finite difference quantity, 
Yh{t,r) in the limit /i — )• admits an asymptotic expansion of the form 

Yh{t, r) = Y{t, r) + h^Cpit, r) + • • • (6.4) 

where h is the discretization scale, ep{t,r) is an /i-independent error function with smooth- 
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ness comparable to Y{t,r), and p is an integer which defines the order of convergence of 
the scheme. Following standard practice, we consider sequences of three calculations per- 
formed with identical initial conditions, but with varying resolutions, h, h/2 and /i/4. We 
then form the differences, ci = Y/i — Y/,/2 C2 = ^/i/2 — ^/i/4) and compute 

log2 (^^^ ~ P- (6.5) 



Fig. 13 shows the results of such a convergence test for the scalar field, ^{t,r), from 
computations in pure harmonic coordinates, and with initial data defined by <I>o — 0.55. 
The plot provides evidence for the expected second order convergence {p = 2) of ^h, and 
similar results are observed for the other dynamical variables. We note, however, that 
there is an obvious degradation of convergence at the highest resolutions used: this issue 
has not been resolved, but may be related to the time-stepping iteration. 



As discussed in Sec. 5.1, in the cases where the Lindblom et al drivers were used to 
evolve the source functions, we used an implicit Euler method to integrate the correspond- 
ing finite difference equations. Since that method is only first-order accurate in time, the 
convergence of the overall scheme in only expected to be first order, and this was in fact 
observed. 



Finally, since we have implemented a free evolution scheme |38|, we can also assess the 
convergence of our numerical solutions by monitoring discrete versions of the Hamiltonian 
and momentum constraints, and M^, respectively. As usual, these constraints are 
defined by contracting the Einstein equations with the unit normal vector to the t = 
const, hypersurfaces, i.e. Ma = n"{Gai3 — T^ii)-, where Gap is the Einstein tensor. In 
order to estimate how well the constraints are satisfied, we discretize them to second 



order, and then compute their L2-norms, as defined by ( |6.2[) , at each time step. Fig. 14 
shows a typical plot of the results for weak initial data (Madm — 0.01) evolved with 
harmonic coordinates. It is clear from the figure that the constraint violations remain 
quite small during the evolution, and that — modulo the previous remark concerning an 
apparent problem at higher resolutions — the constraints are increasingly well satisfied as 
/i ^ 0. 



7. Conclusions 

We have presented a generalized harmonic formulation of the Einstein equations for spher- 
ically symmetric Z)-dimensional spacetimes. Since it is natural to choose coordinates in 
which the symmetries of the geometry are explicit, we have adopted the usual spherical 
coordinates. This results in a coordinate singularity at the origin, r = 0. While at the 
continuum level the equations of motion maintain regularity of a solution which is initially 
smooth at the origin, extra care must be exercised so that this property is reflected in 
discrete numerical calculations. We have thus described a procedure to ensure that the 
origin remains regular in numerical calculations, while preserving the hyperbolicity of the 
evolution system. 
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We have investigated the resulting GH system in the context of fuhy non-Unear gravi- 
tational collapse. To this end we introduced a real, massless scalar field, and have used the 
specification of the initial scalar field profile to control the ensuing strength of the gravi- 
tational interaction. The dynamics that we have considered range from the dispersion of 
weak pulses to the collapse of strong pulses that lead to black hole formation. A key aspect 
of our numerical approach was the use of radial compactification which, in conjunction 
with sufficient dissipation, provided a viable alternative to the truncation of the spatial 
domain and the use of approximate outer boundary conditions. Another ingredient of our 
methodology that was vital for long-term stability of the numerical calculations was the 
addition of constraint-damping terms to the evolution equations. 

Our studies of evolutions using several coordinate drivers lead us to conclude that, in 



spherical geometries, the gauge drivers discussed in |9|, |10|, |18|, |18|] are less effective relative 
to the 3 + 1 simulations that use Cartesian coordinates, and it would be very interesting to 
understand this issue in more detail. Nevertheless, we found that with a certain amount of 
parameter tuning many interesting situations could be successfully simulated with drivers 
that have been proposed in the literature. Perhaps not surprisingly, depending on the 
situation certain drivers performed better than others, leading to longer and/or more ac- 
curate simulations. Specifically, the dynamics of weakly gravitating dispersing pulses could 
be simulated using any of the considered coordinate choices; however the pure harmonic 
gauge arguably provided the cleanest and the simplest choice. For strong-field data, varia- 
tions in the performance of the various drivers were more apparent. In particular, for strong 
but subcritical pulses, the harmonic gauge quickly lead to coordinate pathologies, signaled 
by a collapsing lapse, but this behavior could be partially ameliorated by using one the 
drivers given by (|3.3| ) and (0). The driver (|3.5| ) could also be used to evolve strong- field 
data in some regimes, but the target coordinates which it is designed to asymptotically 
enforce, were not achieved, at least not for the range of the parameters that we explored 
in this work. 

For the case of strong-field, supercritical calculations (i.e. those for which black holes 
form), we found that pure harmonic coordinates could still be of use. In the simulations 
that used excision, it was possible to evolve black holes for as long as a few tens of dynamical 
times. However, the coordinate system remained fairly dynamic even at late times, leading 
to collapse of the lapse near the excision surface on one hand, and to the outwards expansion 
of the coordinate position of the horizon on the other. We were able to use driver conditions 
to moderate the time-dependence, with the best results being obtained through the use of 
the drivers (3.5) with the Bona-Masso target slicing and the F-driven target shift. It would 
be very interesting to find out whether or not parameters and target gauges exist that 
not only slow down the time-dependence of the coordinates at late times, but completely 
eliminate it. 

One of the main goals of this work was to achieve a better understanding of the 
generalized harmonic approach as applied to highly symmetric spacetimes, and to prepare 
ground for an exploration of various gravitational phenomena in axisymmetry using an 
analogous formalism. We expect that the insights gained from our experiments in spherical 
symmetry will also prove useful in the axially symmetric case. In particular, coordinates 
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that are adapted to the axial symmetry are again formally singular on the axis, and the 
equations of motion will need to be regularized there. However, the same regularization 
described above for spherical symmetry can be readily extended to that case. This allows for 
a regular hyperbolic formulation in axial symmetry, which will be discussed in a subsequent 
publication ]3^. 
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A. Asymptotically AdS spacetime 

Here we analyze the asymptotics of AdS spacetime, and discuss a convenient metric ansatz 
as well as a normalization of the source functions. 

The AdSo background can be written in the form, 

ds^ = -(1 + p^/(!^)dT^ + dpV(i + pV^^) + dnl, (A.i) 

where i is the AdS curvature scale. In our model (^) we reproduce asymptotically AdS 
spacetime by letting V{0) = A < that defines f = -{D - 1){D - 2)/A. 

One of the properties of the AdS space is that its asymptotic boundary is time-like: 
in fact, it takes only a finite time for a light signal to propagate to the boundary. Hence, 
in numerical implementations, correct treatment of boundary conditions at spatial infinity 
is crucial. To this end it is useful to transform to conformal coordinates, 

p = £tan{r/£), T = t, (A. 2) 

in which the AdS metric becomes 

ds^ = - cos-^ (r/e) {dt^ - dr^) + f tan^ir/i) d^l- (A.3) 

We note that the entire space has finite extent r G [0,-7r^/2] in these coordinates, but that 
the metric is singular at spatial infinity, r = Tri/2. 

A convenient metric ansatz for evolution using the generalized harmonic approach 
explicitly factors out the background and is given by 

ds"^ = -cos-^ (r/e) gtt {dt^ - dr^) + 2gtrdtdr + f tan^ir/i) e^^dfi^. (A.4) 

In this case the asymptotic behavior of the fields Qab is regular, gab rjab and S" — 0. 
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The source function obtained from (|2.4| ) does not vanish in spherical coordinates even 
in pure AdS where it becomes 

H^"^ = (0, in/p)[l + ((n + 2)/n)(/>V^')]/[l + pVi% (n - 1) cot^i, . . . ,cot0„_i,O) , 

(A.5) 

and where p is given in ( [A. 2 ). In analogy with the asymptotically flat case, we subtract 
a background contribution, which is singular at p = 0, by writing Ha = Ha + Ha'^^ , and 
then use the regular source functions 

/ n 1 + — ''^ \ 
Ha = Htit,r),Hrit,r) + (n - 1) cot ^i, (n - 2) cot 02, • • • , cot . 

V p 1 + F / 

(A.6) 



B. Explicit form of the equations 



We define g2 = gtt drr —Qtr-i to be the determinant of the 2-metric gab-, in (O). The complex 
scalar field is decomposed as <I> = 0^ + ^ • 

The Christoffel symbols and the trace of the extrinsic curvature are given by 

= (gtt g'tr - 9tr g'tt - ^augtr + ^grrgu -ng2Sj/g2 



n / 1 I \ 

Tr = - - + ( -gtt g'rr " '^S'^r g'tt " S'trS'rr + 9rr 5tr - n g2 S' j / g2 



K = a 



ngtr , gtr / I „ Qi ^ ^ ■ I b \ I 
\ grr - 9tr " ^gtr + - grr + n grr S ] g2, 

r ngrr 2 ' 



(B.l) 
(B.2) 



The generalized harmonic equations ( [4.12 ,4.13) in 4L> become 

Rtt — C{t;t) — Ttt = 



7 {9rrf {g'l^ + g'trgrrig'l^ + 5** {gtrf g" - ^g'ld' + g'^g'ttgrrg^' + 



(g'ttY 



^g'^'o'trgW" + ^ + t(5")' {gu f - U 

2g2 4 V 

gtrHt — gtt.Hr gtt\ , , / '^gtt , gttHr — gtrHt 

^ ]9tt+[ — + 

2g2 g2rj 



2 S] + 



52 



grrHt gtrHr gtr \ ■ , / tr\2 I ■ , r, tr tt ! ■ , ^( tr\2 ■ ■ , 

— ]gtt + {g ) gtrgtt + 2g g gugtt + 7;ig ) grrgtt + 



252 



92 r 



g'tr + 
1 



2g''i' gtrgtt -Ht- gtt - ^i'glt - gttV + 2{g"'fg[,gtr, 



(B.3) 
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Rtr — C{t;r) ~ ^tr — 

Ig'^glrg'ttd'^'' - \glrg" + g'^g'^" + \g''g'' + {gtrf - 

^/ . f grrHt - gtrHr gtr\ / , / (fl'*'')^ , rr tt\ I I ^ tr/ , 
gu ^ guHr -gtrHA ^.^ ^ 1 _ ^ 2g^^'gUr + 

g2f / 4 

^ ^.^^.^ ^ + g''-g''g'tr9u + + 



Rrr C{r;r) '^rr 



1 1 29 1 

2(ff")'5ir5« - i^Hr - 2<t>[<i)i - 2ct>',ct>r - iS' S - — - g'^gtr' - ^g^'gir + 

\g''g''grrgtu (B.4) 



Ig'rrg'ttig'''? + HrgrMl^ + + •^g^'grrgU' + s/'-ff;.^;.^*'' + 

g'^g'ttgrrg'^ + 45"5^.5tr<7*^ - + [g'rrf + g"i' (ff*.)' - 

2c?2 \g2r 2g2 J 



grrHt — gtrHr 2gtr\ , , 46" , ( gtr , gtrHr — grrHt \ . , 
]gtr-Hr + + 7y grr + 

52 52?^/ r \g2r 2g2 J 

(g'^fguglr - \ig''r {g'ttf - ^g^'grr - \g''gl. (b.5) 



^66 - C(0;e) — Tee = 

^S'gtt _ _gu_ ^ e"^-^ - V + + H \ + - + 

g2r 52^2 r2 g2r * \^ 52 g2r 52 j 

hJ-^-^ + ^ ) - 2g''S' - g''S - g^'S". (B.6) 

\^ 52 52^ 52 y 

Written in full, the constraint damping terms, Z^^y = k (n(^Cj,) — ^^g^yU^Cp), that we 
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subtract from the above equations to form ( |2.9| ), are 



92 



gtrg'rrQtt^ . 9rrgrrgtt^ gtrHr gtt {grrgtt " '^gtr'^) g'trgtt 



4 52 



+ 



452 



r,/ grrgtrgtrgtt , / 2 

gtr^ gtt 7, ^ [gtr 

2 52 V 

grr {grr gtt — "^gtr"^) gtt 

4^^ 



2 

grrgtt 



+ 



252 

gtr (3 5rr5ft " 45tr^) 5ft 

452 



+ (5rr5ti - 25fr^) S 



(B.7) 



52 



gttgtrgrr"^ _l_ gtrgttgrr^ gttHr 5i 



2 52 



452 



rr _|_ 5ir ^i5 



2 / 

gtt grr grr 



gtrgttg'trgrr . ( grrgtt - '^gtr'^) g'ttgr 



252 



+ 



452 



+ gttS' grr + 



2 452 

gtrgttgrrgrr 



+ 



452 



gtrS grr (B.8) 



gttgtrgrr^ _l_ gtrgttgrr'^ gttHr g 



an 

52 V 2 52 ■ 452 2 
gtrgttg'trgrr , ( 9'rr5it " 25tr^) 5tt5; 



rr ^ 5fr Htg rr 



2 I 

gtt grr grr 



252 



+ 



452 



+ gttS' grr + 



2 452 

gtrgttgrrgrr 



+ 



452 



5tr»S' grr 



(B.9) 



52 



"5tt5r 



4 52 

5tr5tr 5rr 

252 



+ 



Htgrr , gttg'tr9r 



+ 



Sg.r 



2 2 52 

gtrHr 



gtrgttgrr 

452 



5*r 5tt5rr , r" , (^^*'''' 
+ 5tr'J H 



grrgtt gr 



452 



452 



. (B.IO) 



Finally, the Hamiltonian and momentum constraints, Mq = n°'{Gai3 — T^/s), take the 
form 



Mt = -^g2gtt {(t>'i)^ - ^g2gtt {(prY - ^g2gtt {S'Y + ^525rr {4>i) + ^525rr [(Vr 



,„,2 52(52^2-6 2'^52+5ti) „ c.// , 4 (5rr5tr5it - 5ir^) 5' 
525rr [S] ^ 2 ~ ^92gttS H ^ ^ 



/ f rr2o' 2 I rr tr a 2 , gtt^\ <-,/ / ^g2gtt o / 3 , 2 rr tr\ c'^ 1 
5rr ( 5 5^52 + 5 5 ^52 + — 1 + 5 ( 8 (5t/ + 52^5rr5 5 ) I + 

/ fn rr tr o' 2 , r, rr tt 2 '^gtrgtt\ 1 ( tr"^ o' 2 , tr tt a 2 , gtr^\ 

gtri^g g 5" 52 +25 5 Sg2 — I +5**1 5 5^52 +5 5 '^52 +~) + 

{2grrgtrgtt - '2gtr^) S' - g2grrS. (B.ll) 
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Mr = [gtr^ + g2^9rrg'''g'') (00' + '^929rM[ + {gtr' + g2^grr9"g'') (Kf 

^g2gtrS" + (-S'g2 - —) g'rr + 2g2grr4'r'i'r + '^^'^^"'^ + 2g2grrS' + 



r,/ / '^{gtr^ - grrgtrgtt) c\ ^ ' ( rr tr qI 2 tr'^ d 2 , 9tr9tt\ . 
S — + 'ig2grrS + ffrr -g 9 S 52 " 9 892^ H + 



9[^ ^-2g'^'s'g^ - Ig'^g^'Sg^^ - ^) - 2^25*. [S' f + 

9'u {-9''9''S'g2^ - i''Sg2' + ^) . (B.12) 
C. Discretization 

The second order accurate finite difference approximations (FDAs) for the time derivatives 
on a uniform grid with spacings Ar, At at a point (n, i) (see Fig. [l[) are 

- 2A^ ' 

yn+l _ 2 yn , y"--! 

a^y." = ^ \/ ^ — . (c.i) 



Here "second order" means that the continuum expression is approached by the FDA 
counterpart at a rate O(At^). For the spatial and mixed derivatives the stencil is modified 
depending on the position of the mesh point relative to the extremities of the grid. We use 
second order accurate expressions of the form 

• Centered derivative. 

yn _ yn 

~ 2Ar ' 



(Ar)2 

-] 

4 Ar At 



yn+l _ yn— 1 _ yn+l , yn— 1 
52^yn ^ ^^+1 +^^-1 (C.2) 



One-sided (backward) derivative. 

' ~ 2Ar 

2..n ^yi'-5y:ii+^Yr+2-y:i3 



(Ar)2 
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